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Abstract 

Background: Flux coupling analysis (FCA) has become a useful tool in the constraint-based analysis of genome-scale 
metabolic networks. FCA allows detecting dependencies between reaction fluxes of metabolic networks at 
steady-state. On the one hand, this can help in the curation of reconstructed metabolic networks by verifying whether 
the coupling between reactions is in agreement with the experimental findings. On the other hand, FCA can aid in 
defining intervention strategies to knockout target reactions. 

Results: We present a new method F2C2 for FCA, which is orders of magnitude faster than previous approaches. As a 
consequence, FCA of genome-scale metabolic networks can now be performed in a routine manner. 

Conclusions: We propose F2C2 as a fast tool for the computation of flux coupling in genome-scale metabolic 
networks. F2C2 is freely available for non-commercial use at https://sourceforge.net/projects/f2c2/files/. 



Background 

The huge amount of genomic, transcriptomic and related 
data has allowed for a fast reconstruction of an increasing 
number of genome-scale metabolic networks, e.g. [1-7]. 
In the absence of detailed kinetic information, constraint- 
based modeling and analysis has recently attracted 
ample interest due to its ability to analyze genome-scale 
metabolic networks using very few information [8-10]. 
Constraint-based analysis is based on the application of 
a series of constraints that govern the operation of a 
metabolic network at steady state. This includes the sto- 
ichiometric and thermodynamic constraints, which limit 
the range of possible behaviors of the metabolic network, 
corresponding to different metabolic phenotypes. Apply- 
ing these constraints leads to the definition of the solution 
space, called the steady-state flux cone [11]: 



C = {v E ' 



Sv = 0, Vi > 0, for all i e Irr}, 



(1) 



where S is the m x n stoichiometric matrix of the net- 
work, with m internal metabolites (rows) and n reactions 
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(columns), and the vector veR" gives a flux distribution. 
Furthermore, Irr c {1, ...,n} denotes the set of irre- 
versible reactions in the network, and Rev = {1, . . . , n}\Irr 
denotes the set of reversible reactions. 

The flux cone contains the full range of achievable 
behaviors of the metabolic network at steady state. Var- 
ious approaches have been proposed either to search 
for single optimal behaviors using optimization-based 
methods [12-16] or to assess the whole capabilities of a 
metabolic network by means of network-based pathway 
analysis [11,17-20]. 

Flux coupling analysis (FCA) is concerned with describ- 
ing dependencies between reactions [21]. The stoi- 
chiometric and thermodynamic constraints not only 
determine all possible steady-state flux distributions over 
a network, they also induce coupling relations between 
the reactions. For instance, some reactions may be unable 
to carry flux under steady-state conditions. If a non-zero 
flux through a reaction in steady-state implies a non-zero 
flux through another reaction, then the two reactions are 
said to be coupled (see Def. 2 for a formal definition). FCA 
has been used for exploring various biological questions 
such as network evolution [22-24], gene essentiality [22], 
gene regulation [25-27], analysis of experimentally mea- 
sured fluxes [28,29], or implications of the structure of 
the human metabolic network for disease co-occurrences 
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[30]. Having a time efficient implementation of FCA is 
important in such studies. 

After introducing the main existing algorithms for flux 
coupling analysis, we propose in this paper a new algo- 
rithm which significantly speeds up the calculation of 
flux coupling. Our algorithm is based on two main 
principles. First, we reduce the stoichiometric model 
as much as possible when parsing the stoichiomet- 
ric matrix. Second, we use inference rules to mini- 
mize the number of linear programming problems that 
have to be solved. We prove the efficiency of our algo- 
rithm by successfully competing with the most recent 
approach [31]. We show that FCA can now be quickly 
performed even for very large genome-scale metabolic 
networks. 

Approaches for flux coupling analysis 

Several algorithms were developed to calculate flux cou- 
pling between reactions. For a comparison among the 
existing approaches, the reader may refer to [31,32]. In the 
following, we focus on flux coupling methods based on 
solving a sequence of linear programming (LP) problems. 
These methods have proved to be significantly faster than 
other algorithms. 

Definitions 

We give a short overview of the important concepts we 

will use throughout this paper. First, we formally define 

blocked reactions in a metabolic network. 

Definition 1 (Blocked reaction). Given the steady-state 

flux cone C, let i e {1, ...,n} be a reaction. If Vi = 

0, for all v e C, reaction i is called blocked, otherwise i is 

unblocked. 

In the following, we assume that the flux cone is not triv- 
ial, i.e., not all reactions are blocked. Next, we define the 
(un)coupling relationships between reactions. 
Definition 2 (Coupling relations). Let if be two 

unblocked reactions. The (uncoupling relationships 

— <>, ^ and i-> are defined in the following way: 

• i ^ j if for all v e C, v; = 0 implies vj = 0. 

• i j if for all v e C, v; = 0 is equivalent to vj = 0. 

• i ^ j 

if there exists X ^ 0 such that for all v e C, vj = Xvi. 

• i i-> j if there exists v e C such thatvt = 0 and vj ^ 0. 

Reactions i and j are fully (resp. partially, directionally) 

coupled if the relation i ^ j (resp. i j, i -> j) holds. 
Otherwise, i and j are uncoupled. 

Note that i ^ ; (resp. i <+ j) is equivalent to j ^ i (resp. 

j i). In addition, i ^ ; implies i which in turn is 

equivalent to (i — £ j and ; — £ i). 



As shown in [33] the reversibility type of reactions is a 
key concept in flux coupling analysis. 
Definition 3 (Reversibility types). A reversible reaction 
i e Rev is called fully reversible if there exists a flux vector 
v e C such that V{ ^ 0 and vj = 0 for all j e Irr. Otherwise, 
reaction i is called pseudo-irreversible. 

Using the reversibility type of reactions, we define the 
following reaction sets which will be used to determine 
the cases where flux coupling can occur between reac- 
tions. 

• Frev = {i | i is fully reversible}, 

• Prev = {i | / is pseudo-irreversible and there exist 
v + , v" e C such that v+ > 0, vj" < 0}, 

• Irev = {i | i £ Frev U Prev and v; ^ 0 for some v e 
C}, 

• Blk = {i | / is blocked}. 

Note that the above reaction sets are disjoint and their 
union is equal to the set of all reactions, i.e., Blk U Irev U 
Prev U Frev = {1, . . . , n}. 

As the classification of reactions according to their 
reversibility types has been treated elsewhere [31,33,34], 
we focus on calculating the flux coupling between reac- 
tions given the reaction sets Blk, Irev, Prev and Frev. First, 
we briefly review the main existing flux coupling methods 
based on linear programming. Afterwards, we present our 
new method to speed up the flux coupling analysis. 

Flux coupling finder 

The Flux Coupling Finder (FCF) algorithm [21] deter- 
mines blocked reactions as well as coupled reactions by 
solving a sequence of linear programming (LP) problems. 
The FCF algorithm requires that each reversible reac- 
tion is split into two irreversible reactions, one forward 
and one backward. This may hamper the application of 
FCF to determine flux coupling in large genome-scale 
metabolic networks. Indeed, splitting reversible reactions 
results in an increase in the number of variables (resp. 
constraints) by \Rev\ (resp. 2\Rev\). Since the FCF algo- 
rithm solves four LP problems to identify the maximum 
and minimum flux ratios for every pair of reactions, the 
total number of LP problems that have to be solved 
is very large. Furthermore, the FCF algorithm does not 
compute directly coupling relationships between reac- 
tions. A post-processing step is needed to deduce cou- 
plings between reactions (in the original network) from 
those between reaction directions (in the reconfigured 
network). More importantly, the FCF algorithm explores 
exhaustively all possible reaction pairs. This leads to a very 
big number of LP problems that have to be solved. This 
strategy may not scale well for genome-scale models of 
complex microorganisms which involve a large number 
of reactions. 



Larhlimi etal. BMC Bioinformatics 2012, 13:57 
http://www.biomedcentral.eom/1 471 -21 05/1 3/57 



Page 3 of 9 



Reversibility-based flux coupling analysis 

To cope with the drawbacks of the FCF algorithm, we used 
the reversibility type of reactions to develop an improved 
version (WRP-FCF) of the original FCF method [31,34]. 
When looking for coupled reactions, WRP-FCF applies 
linear programming only in those cases where coupling 
relationships can occur [33]. Namely, given two unblocked 
reactions i and we have exactly three cases where a flux 
coupling is possible between i and 

Flux coupling between reversible reactions 

If i,j e Prev or i,j e Frev, i ^ j,i j and i ^ ; 
are equivalent. More importantly, only the stoichiometric 

constraints determine whether i ^ j holds, indepen- 
dently of the thermodynamic constraints. 

Assume that blocked reactions have been identified 
beforehand and their respective columns in the stoichio- 
metric matrix have been removed. Consider the following 
LP problem 

max{vj : Sv = 0, v { = 0, 0 < vj < 1}, (2) 
and let v* be an optimal solution. According to Proposi- 
tion 6.13 in [34], i ; if and only if v* = 0. 

Flux coupling between (pseudo-)irreversible reactions 

If i e Irev and ; e Prev, we only have to check whether 

i ^ The other coupling relationships cannot occur. Let 
v 1 resp. v 2 be optimal solutions of the two LP problems 

max{v/ : Sv = 0, v^ > 0, k e Irr, v/ = 0, Vj < 1}, 
min{v/ : Sv = 0, v^ > 0, k e Irr, V{ = 0, Vj > — 1}. 

(3) 

Then i -> j if and only if vj = vj = 0. 

Flux coupling between irreversible reactions 

If g Irev, in analogy with the FCF algorithm, we 
determine upper and lower bounds Uij and such that 
0 < LijVj < Vi < UijVj for all v e C by solving the LP 
problems 

Uj = min {v/: Sv = 0, vj = 1, v k > 0, k e Irr}, 
Uij = max [vf. Sv = 0, Vj = 1, v/ c > 0, k e Irr}. 

Comparison of Ly and Uij allows us to determine 
whether reactions i and ; are coupled using the following 
rules: 

• i ^ ; (resp. ; ^ /) if and only if Lg i=- 0 (resp. 
Uij £ +oo), 

• j ^ i if and only if Ly ^ 0, Uij ^ 0 and Ly = Uij. 

The improved version WRP-FCF does not make an 
exhaustive computation for all pairs of reactions and does 
not require a reconfiguration of the metabolic network. 



Accordingly, not only is the number of LP problems used 
by WRP-FCF smaller than the number solved by FCF, 
but also the LP problems used by WRP-FCF are simpler 
than those employed by FCF. For mathematical proofs, the 
reader may refer to [31,34]. 

Feasibility-based flux coupling analysis 

Linear programming can be solved by enumerating a set 
of adjacent feasible solutions such that at each solution the 
objective function improves or is unchanged. Accordingly, 
searching for a feasible solution can be faster than com- 
puting an optimal solution [35]. Based on this observation, 
[31] proposed the FFCA approach for flux coupling anal- 
ysis transforming the optimization-based LP problems 
used in the WRP-FCF method into feasibility-based LP 
problems. The authors compared FFCA with other avail- 
able flux coupling algorithms, and showed that FFCA is 
slightly faster than WRP-FCF. 

Methods 

Before using linear programming to calculate flux cou- 
pling between reactions, we preprocess the metabolic 




(a) Exemplary metabolic network MetNet 



'System boundary" | 8,9 




(b) MetNet after preprocessing 



Figure 1 Exemplary metabolic network MetNet before and after 
preprocessing, (a) MetNet consists of eight metabolites (A, ...,H), 
and thirteen reactions (1 , . . . , 1 3), whereof six reactions are 
irreversible. Metabolites are depicted as nodes while reactions are 
depicted as arrows. Reversible reactions are indicated by double 
arrowheads, (b) MetNet after preprocessing where dead-end 
metabolites and blocked reactions were removed and fully coupled 
reactions were merged iteratively. This resulted in the removal of the 
blocked reaction 13 and the merging of the pairs of equivalent 
reactions (1,2), (8, 9) and(1 1, 12)). 
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Table 1 Main steps of the preprocessing procedure 

Step Rule 

1 . Identify dead-end metabolites and the corresponding blocked reactions. 

2. Apply Reduction Rule 1 to remove the rows (resp. columns) corresponding to dead-end metabolites (resp. blocked reactions) from the 
stoichiometric matrix. 

3. Apply the TFC rule to determine reactions which are proportional to each other and update their reversibility type. 

4. Apply Reduction Rule 2 to keep only one column for each set of reactions which are proportional to each other. 

5. Repeat Steps 1 -4 until neither a dead-end metabolite nor a pair of fully coupled reactions can be identified. 

6. Update the reversibility type of each reaction and remove the columns corresponding to blocked reactions from the stoichiometric 
matrix. 

7. Use a basis of the kernel of the stoichiometric matrix to identify fully coupled reactions. This step is optional as the kernel computation 
may lead to numerical problems. 

8. Classify reactions according to their reversibility type. 



network in order to i) reduce the number of variables and 
constraints of the subsequent LP problems and ii) classify 
reactions according to their reversibility type. The net- 
work reduction is mainly based on the removal of trivially 
blocked reactions and the merging of the stoichiomet- 
ric columns corresponding to trivially coupled reactions 
[36-39]. For this, one can use the kernel of the stoichio- 
metric matrix. Alternatively, we can apply the following 
reduction rules which require a simple parsing of the 
stoichiometric matrix and are not time consuming. This 
strategy allows avoiding potential numerical instabilities 
related to the computation of a basis of the kernel. 

Preprocessing 

Certain metabolites, called dead-end metabolites [37], are 
produced (resp. consumed) by some reactions without 
being consumed (resp. produced) by other reactions. This 
concept is illustrated in Figure 1(a) where the dead-end 
metabolite H is produced by reaction 13 without being 
consumed by any of the remaining reactions. 

As stated below, reactions which are consuming or pro- 
ducing dead-end metabolites are blocked. 
Observation 1 (Dead-end Metabolite). Letk e {1, . . . , m} 
be a metabolite. Then, the following hold: 

• If there exists a reaction i such that S# 7^ 0 and 
Sty = 0 for all j ^ i, then reaction i is blocked. 



Table 2 Transitivity inferred flux (un)coupling 



Known flux coupling 


i ^ / 


. =o . 

1 J 


. =o . 
i^j 


. =o . 
/ ~+ i 


k ^ / 


k^j 


. =0 . 


k^j 


j^k 




k o / 


k 




k^j 




k -> / 


k^j 


k^j 


k^j 




k \-> i 


k \-> j 


k\-> j 




k \-> j 


i \-> k 


j \-> k 


j \-> k 


j \-> k 





• If there exists a set of reactions I c Irr such that 
Ski > 0 (resp. S^i < 0) f° r all i £ I and Sq = 0 for all 
j £ I, then all reactions i e I are blocked. 

In each of these cases, k is called a dead-end metabolite. 

Certain metabolites are involved in exactly two reac- 
tions. For instance, in the network MetNet depicted in 
Figure 1(a), metabolite E is produced/consumed only by 
reactions 8 and 9. The following observation states that 
the fluxes through reactions involving such metabolites 
are proportional to each other. 

Observation 2 (Trivial Full Coupling (TFC)). Let i and 
j be two reactions such that, for some metabolite k e 
{!,..., m}, Sj^ 7^ ^ Sjg 7^ 0 and S^i = 0 for all I ^ i,j. 
Then, reactions i and j are either blocked or fully coupled. 

The identification of dead-end metabolites and their 
corresponding blocked reactions allows us to reduce the 
number of metabolites and reactions that matter for 
identifying coupled reactions. As stated in the follow- 
ing observation, the removal of the rows (resp. columns) 
in the stoichiometric matrix corresponding to dead-end 
metabolites (resp. blocked reactions) does not influence 
the flux coupling between reactions. 
Observation 3 (Reduction Rule 1). Let D be a set of dead- 
end metabolites and let B be a set of blocked reactions. For 
convenience, suppose B = {s + 1, . . . , n}. Let S r be the sub- 
matrix of S formed by the rows with k £ D and the 
columns S*i with I £ B. Let Irr' = Irr\B. Then, for all pairs 
of reactions i £ B and j £ B, 

• i ^ / if and only if v\ = 0 implies v'- — 0, for all 

v' e W such thatS'v' = 0andV p >0 for all p e Irr'. 

• i ^ j if and only if there exists A/ ^ 0 such that 
v'. = k V., for all v' e W with SV = 0 and 

v' p > 0 for all p e Irr' . 

The next observation shows that two fully coupled reac- 
tions can be represented by only one column in the 
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Table 3 The main steps of the F2C2 algorithm 

Step Rule 

1 . Apply the preprocessing procedure shown in Table 1 . 

2. Apply the feasibility rule using the feasible solutions obtained when solving the LP problems used in Step 1 . 

3. Apply the TDC and TUC rules to determine trivially (un)coupled reactions. 

4. Identify fully coupled reversible reactions by solving the LP problems (2). This is not necessary if the kernel of the stoichiometric matrix is 
used in Step 1 . 

5. Determine the dimension t of kern(S*R ev ) and remove t independent fully reversible reactions and their coupled reactions. This step is 
optional since t is often small. 

6. Determine the flux coupling between (pseudo-)irreversible reactions by solving the LP problems (3) and (5). 

7. For each LP problem solved in Step 6, use the inference rules in Table 2 in combination with the feasibility rule. 



stoichiometric matrix, without altering the flux coupling 
between reactions. 

Observation 4 (Reduction Rule 2). Let k, I be two reac- 
tions such that for all v e C, v/ = Xv^ for some X ^ 0. 
For convenience, suppose I = n and X > 0. Let S r be the 
mx(n — l) matrix defined by S^ p = S* p for all p ^ k,l and 
S f ^ k = S* k + A.S*/. Let In* = (Irr U {/<}) \ {/} if I e Irr, and 
Irr' — Irr otherwise. Then, for all pairs of reactions i ^ I 
and j 7^ l } 

• i — £ j if and only if i/ = 0 implies v'- = 0, for all 
v' e R"- 1 such thatS'v' = 0 and 

v p >0 for all p e Irr' . 

• i ^ j if and only if there exists X f ^ 0 such that 

v'. = A. V., for all v' e R"' 1 with S f v f = 0 and 
j 1 

v' p > 0 for all p e Irr' . 

Note that when applying the reduction rules of Obser- 
vations 3 and 4, further metabolites and reactions may 
fulfill the conditions of Observations 1 and 2. Accord- 
ingly, we apply these reduction rules in an iterative way. 
As an illustration, the reduction of the network Met- 
Net depicted in Figure 1(a) involves two iterations. In 
the first one, metabolite H and reaction 13 are removed, 
the pairs of reactions (1,2) and (8,9) are merged and 
metabolites A and E are removed. In the second iteration, 
the equivalent reactions (11, 12) are merged and metabo- 
lite G is removed. The preprocessed network depicted 
in Figure 1(b) contains only four metabolites and nine 
reactions. 

Certain fully coupled reactions could not be identified 
using Observation 2. The following lemma proves that all 
fully coupled reaction pairs can be deduced from the ker- 
nel kern(S) = {v e R n \ Sv = 0} of the stoichiometric 
matrix after the removal of all blocked reactions. 
Lemma 1. Let (S,Irr) be a metabolic network with n 
unblocked reactions. For a pair of reactions (/,;) the follow- 
ing are equivalent: 

• i ^ /. 



• there exists X e R \ {0} such that v,- = Xvj, for all 
v e kern(S). 

Proof. <^=" Immediate. 

Since i ^ there exists X ^ 0 such that v; = Xvj 
for all v e C. Assume by contradiction that there is v e 
kern(S) such that v; ^ Xvj and let L = {I e Irr | v\ < 0}. 
Since every reaction is unblocked, for every / e L there 
exists^ e C with = 1. Let w = v — J2ieL v l8^- 
Clearly, w e kern(S) and w\ > 0 for all / e Irr, thus w e 
C. However, W{ ^ Xwj, contradicting i ^ j. The required 
statement follows. □ 

In analogy with the WRP-FCF and FFCA approaches, 
we identify the reversibility type of reactions in order 
to apply linear programming only in those cases where 
coupling relationships can occur [33]. Here, we use the 
procedure for reaction classification described in [31,34]. 
Note that applying the above reduction rules beforehand 
reduces the number of variables and constraints in the LP 
problems used for the classification of reactions. 

Based on the results above, we propose to apply the pre- 
processing procedure given in Table 1 before identifying 
coupled reactions using linear programming. We show 



Table 4 Genome-scale metabolic networks with the 
number of metabolites (tfmet.) and reactions (tfreac.) 
before and after applying the preprocessing steps 





Original size 


Prepr. size 


Network name 


tfmet. 


tfreac. 


ttmet. 


ttreac. 


M. barkeri, iAF692 


628 


690 


149 


221 


S. cerevisiae, /7VD750 


1061 


1266 


248 


446 


M. tuberculosis, /7VJ661 


826 


1025 


240 


418 


E. coll, UR904 


761 


1075 


269 


565 


E. coll, iAF] 260 


1668 


2382 


604 


1272 


E. coll, UO] 366 


1805 


2582 


651 


1376 


H. sapiens, Recon] 


2766 


3742 


725 


1573 
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Table 5 Performance comparison between the FFCA and F2C2 algorithms in terms of the number of LP problems solved 
(#LPs) and their total running times (TRT) 







FFCA 




F2C2 




Nptwnrk 

INCl VV \J 1 l\ 




TRT 




TRT 


Prepr. RT 


M. barken, iAF692 


301975 


59 m 40 s 


774 


7s 


5s 


S. cerevisiae, /'/VD750 


472629 


1 h 50 m 17s 


1280 


21 s 


15s 


M. tuberculosis, /YVJ661 


556504 


3 h 5 m 36 s 


1506 


22 s 


16s 


E. coli, UR904 


655437 


2 h 40 m 33 s 


1580 


26 s 


18s 


E. coll, iAF] 260 


4256786 


4 d 31 m 26 s 


3309 


2 m47s 


2m 


E. coli, UO] 36. 


4877262 


4 d 5 h 30 m 46 s 


3955 


3 m 55 s 


2 m 45 s 


H. sapiens, Reconl 


4566304 


4d 18h3m37s 


3903 


5 m20s 


4m9s 



For the F2C2 algorithm, TRT includes the time (Prepr. RT) required for the preprocessing step. Computation times are given in days (d), hours (h), minutes (m) and 
seconds (s). 



later that the preprocessing step turns out to be crucial for 
obtaining an efficient flux coupling algorithm. 

Algorithmic improvements 

In certain metabolic networks, the conversion of a set 
of substrates into a set of products can be made by dif- 
ferent reactions having the same stoichiometry. A sim- 
ple example of such reactions are isoenzymes which 
make the same conversion of substrates into prod- 
ucts. This concept is illustrated in Figure 1(a) where 
both reactions 4 and 5 convert C into D in the same 
way. This holds also for reaction 7 and the merged 
equivalent reactions (8,9) in Figure 1(b) showing that 
the network preprocessing may simplify the identifi- 
cation of such alternative conversions. The flux cou- 
pling of such reactions is trivial using the following 
lemma. 

Lemma 2 (Trivial Uncoupling (TUC)). Let i,j e Irev and 
k, I e Prev U Frev be four reactions. Then, the following 
holds: 

• If S*i = aS*j for some a > 0, then i i-> p and j i-> p 
forallp^Elk. 

• If S*i = aS*j for some a < 0, then p h-> i (resp. 
p h> j) for all p <£ Elk U {/} (resp. p <£ Elk U {/} ). 

• If S*i = otS*k for some a^0, then i h-> p and p h-> i 
for all p £ Elk and q \-> k for all q £ Elk U {/}. 

• If S*k = aS*i for some a ^ 0, then /ck p and p \-> k 
for all p £ Elk U {/} and I \-> q and q \-> I for all 

q <£ Elk U {/<}. 

Proof. The proofs of the four statements are similar. We 
only consider the first one. Suppose S*i = aS*j for some 
a > 0 and let us prove i \-> p. Let p £ Elk. There exists 
v G C such that v p ^ 0. Let w e W such that w t = 0, 
Wj = avi + vj and w q = v q for all q ^ i,j. We have w e C, 
W( = 0, w p 7^ 0 and so the claim follows. □ 



The next observation states that metabolites which are 
involved only in irreversible reactions and consumed or 
produced by exactly one reaction define trivial directional 
couplings between these reactions. 

Observation 5 (Trivial Directional Coupling (TDC)). Let 
k be some metabolite such that Sj^ = 0 for all I e Frev U 
Prev. Let P = {i | S ki > 0} and N = {/ | S kj < 0}. If 

card(P) = 1 (resp. card(N) = 1), then i —> j (resp. j — £ i) 
for all (/,;) e P x N. 

Since directional flux coupling is a transitive relation, 
the flux (un) coupling between many reactions can simply 
be deduced from dependencies between reactions whose 
flux (un)coupling has been determined beforehand. This 
allows us to significantly reduce the total number of LP 
problems to be solved. Examples of such inferred flux 
(un)couplings are given in Figure 1(b). According to the 

TDC rule, we have (1, 2) (8, 9). By solving the LP prob- 
lems (4), we obtain 10 \-> (8,9). We can easily conclude 
that 10 h> (1,2). 

Table 2 shows the inferred flux (un) coupling relations 
we can deduce from known relationships between reac- 
tions. 

For some pairs of reactions, we need to solve at least one 
LP problem. The optimal solution not only determines the 
flux coupling between the considered pair of reactions, 
but also allows one to determine many other uncoupled 
reactions. 

Observation 6 (Feasibility Rule). Let v e C be a steady 
state flux vector and let I = {i \ V[ — 0} and J = {j \ Vj ■ / 
0}. Then i \-> j for all (/,;) e I x /. 

In general, most irreversible reactions are uncoupled 
to each other. Accordingly, the LP problems (4) used 
to determine coupled irreversible reactions are often 
unbounded. This limits the use of the feasibility rule, 
which requires the calculation of a feasible flux vector. In 
order to optimally use the feasibility rule, instead of solv- 
ing the LP problems (4) to decide whether two irreversible 
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Figure 2 Visualization of the LP problems solved using different 
algorithms, (a) The FFCA algorithm, (b) the FFCA algorithm after 
applying the preprocessing procedure given in Table 1 without 
kernel computation (Step 7), (c) the FFCA algorithm after applying the 
preprocessing procedure and using the kernel of the stoichiometric 
matrix to identify fully coupled reactions and (d) the F2C2 algorithm 
given in Table 3. The dashed lines mark the boundary of the I rev, Prev 
and Frev regions. Colors: Black - an LP problem is solved for the 
corresponding (ordered) pair of reactions; Gray - the corresponding 
LP problem is not solved due to one (or both) reactions being 
eliminated from the network (a preprocessing improvement); White - 
corresponding LP problem does not get solved due to an algorithmic 
improvement. 



reactions e Irev are coupled, we propose to solve the 
bounded LP problems 

Lij = min{v; : Sv = 0, Vj = 1, > 0, k e Irr}, , . 
Lji = min{v/ : Sv = 0, v; = 1, > 0, k e Irr}, 

and to use the following scheme: 

• i ^ / (resp.; ^ i) if and only if ^ 0 (resp. Lji ^ 0), 

• j ^ i if and only if Lij ^ 0, Lji ^ 0 and Lij = I /Lji. 

The following observation states that removing a fully 
reversible reaction does not alter the flux coupling 
between (pseudo-)irreversible reactions. 
Observation 7. Let k e Frev be a fully reversible reaction. 
For convenience, suppose k = n. Let S f be the m x (n — 1) 
matrix defined by S^ p = S* p for all p ^ k and let Irr' = Irr. 
Then, for all pairs of reactions i £ Frev and j £ Frev, 

• i ^ j if and only if v\ = 0 implies v'- = 0, for all 

v' e W- 1 with S f v f = 0andv f p >0 for all p e Irr'. 

• i ^ j if and only if there exists A/ ^ 0 such that 
vj = A.V., for all v' e IT" 1 with S f v' = 0 and 
v' p > 0 for all p e Irr' . 

Let S*R ev be the submatrix denned by the columns in 
S corresponding to the reversible reactions and let t be 
the dimension of kern(S*R ev ). Based on Observation 7, we 
can remove up to t independent fully reversible reactions 
without altering the flux coupling between (pseudo-)- 
irreversible reactions. Since certain fully reversible reac- 
tions may change their reversibility type after the deletion 
of a fully reversible reaction, we first remove a randomly 
chosen reaction i e Frev together with the coupled reac- 
tions with L We calculate the impact of this deletion on 
the dimension of kern(S*R ev ). If this dimension decreases 
by 1, the deletion is maintained; otherwise the removed 
reactions are restored to the network. This is repeated 
until t independent fully reversible reactions and their 
coupled reactions are removed. We assume that the flux 
coupling between fully reversible reactions is determined 
beforehand. 
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Based on the above results, we propose the Fast Flux 
Coupling Calculator (F2C2) to determine coupled reac- 
tions. The main steps of the F2C2 algorithm are given in 
Table 3. 

Results and discussion 

The F2C2 algorithm has been implemented within the 
MATLAB environment, using CLP (via the Mexclp 
[40] interface) as the underlying linear programming 
solver. For benchmarking, we analyzed the following 
genome-scale metabolic networks: Escherichia coli, 
iJR904 [1], Saccharomyces cerevisiae, iND750 [2], 
Methanosarcina barkeri, iAF692 [3], Mycobacterium 
tuberculosis, /A//661 [4], Escherichia coli, /AF1260 [5], 
Homo sapiens, Reconl [6] and Escherichia coli, //01366 
[7]. For the numerically sensitive parts, a tolerance level of 
1 0 e - 6 was set. All computations were performed using 
a single Intel Xeon 5160 (3.0 GHz) processor on a 64-bit 
Debian Linux 6.0 system. 

As pointed out in the previous section, part of the per- 
formance gain of F2C2 over previous FCA algorithms 
stems from the fact that the preprocessing steps reduce 
the network size. This affects the running time on two lev- 
els: there are fewer reaction pairs and the LP problems to 
be solved have reduced size. The dramatic effect of the 
preprocessing steps on the network size is presented in 
Table 4. 

The algorithmic improvements further reduce the num- 
ber of LP problems to be solved. A direct performance 
comparison between the FFCA and F2C2 algorithms 
(including the running times and number of LP prob- 
lems solved) is summarized in Table 5. In all cases, F2C2 
outperformed FFCA by several orders of magnitude. In 
[31] it has been shown that FFCA is more efficient on 
genome-scale metabolic networks than other flux cou- 
pling algorithms. 

For an intuitive, visual representation of the individual 
improvements' impact on the algorithms performance, 
a more in-depth analysis has been performed on the 
recent metabolic network of E. coli, //01366. Four differ- 
ent sets of improvements were cumulatively switched on, 
and the linear programs solved were plotted for each case 
(Figure 2). To better highlight the relevant differences, the 
following modifications were applied to the results. First, 
249 (out of 2582) reactions identified as blocked were 
removed from the images. This is a common preprocess- 
ing step in most FCA algorithms. Secondly, the order of 
reactions was permuted such that the reactions in Irev, 
Prev and Frev are clustered together. Additionally, in each 
of these three clusters, the fully coupled reactions were 
moved towards the end of the segment. 

Figure 2(a) plots the LP problems solved in the FFCA 
algorithm. Applying the simple preprocessing steps with- 
out using the kernel (Figure 2(b)), several reactions are 



found to be fully coupled with others, and as such can be 
merged together. When Lemma 1 is applied (Figure 2(c)), 
all fully coupled sets are detected. As a consequence 
the gray stripes get thicker and the LP problems corre- 
sponding to (Prev, Prev) and (Frev, Frev) pairs need not be 
solved anymore. The use of the algorithmic improvements 
(Figure 2(d)) filters the pairs in (Irev, Irev) and (Irev, Prev) 
blocks, greatly reducing the total number of LP problems 
solved. 

Conclusions 

We have presented the new flux coupling algorithm F2C2, 
which outperforms previous methods by orders of mag- 
nitude. Flux coupling analysis of genome-scale metabolic 
networks can now be performed in a routine manner. A 
software tool F2C2 is freely available for non-commercial 
use at https://sourceforge.net/projects/f2c2/files/. 
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